Namn och CID på gruppmedlemmar: Blend Ahmed Omar (blend), Albin Östling (ostlinga), Ebbe Ledin (ebbel)
import numpy as np # Standard paket för att hantera matamatik och arrayer
import matplotlib.pyplot as plt # Standard paket för att plotta figurer
plt.style.use("ggplot")
import scipy.io as io # Vi lånar funktionen 'io' från scipy för att smidigt kunna importera .mat-filer
# Funktioner för HUPP:en
def fft2c(x):
'''
2D Fourier transform
Denna är perfekt som den är. Bara att använda!
'''
return np.fft.fftshift(np.fft.fft2(np.fft.fftshift(x)))
def ifft2c(x):
'''
2D inverse Fourier transform
Denna är perfekt som den är. Bara att använda!
'''
return np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(x)))
def PAS(E1, L, N, a, lam0, n):
'''
Funktion för att propagera E1 sträckan L genom PAS
'''
# Varje sampelpunkt i k-planet motsvarar en plan våg med en viss riktning [kx,ky,kz]
delta_k =2*np.pi/(N*a) # Samplingsavstånd i k-planet
kx = np.arange(-(N/2)*delta_k, (N/2)*delta_k, delta_k) # Vektor med samplingspunkter i kx-led
ky = kx # och ky-led
KX, KY = np.meshgrid(kx,ky) # k-vektorns x- resp y-komponent i varje
# sampelpunkt i k-planet
k =2*np.pi*n/lam0 # k-vektorns längd (skalär) för en plan våg i ett material med brytningsindex n
KZ = np.sqrt(k**2-KX**2-KY**2, dtype=complex) # k-vektorns z-komponent i varje sampelpunkt.
# dtype=complex tillåter np.sqrt att evaluera till ett komplext tal
fasfaktor_propagation = np.exp(1j*KZ*L) # Faktor för varje sampelpunkt i k-planet
# multas med för att propagera sträckan L i z-led
A = a**2/(2*np.pi)**2*fft2c(E1) # Planvågsspektrum i Plan 1
B = A*fasfaktor_propagation # Planvågsspektrum i Plan 2 (Planvågsspektrum i Plan 1 multat med fasfaktorn för propagation)
E2 = N**2*delta_k**2*ifft2c(B)
return E2
# Propagera fält med PAS
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'
N = 2**10 # NxN är antalet samplade punkter (rekommenderad storlek N=1024)
sidlaengd_Plan1 = 4e-3 # Det samplade områdets storlek (i x- eller y-led) i Plan 1 (rekommenderad storlek 4 mm)
a = sidlaengd_Plan1/N # Samplingsavstånd i Plan 1 (och Plan 2 eftersom vi använder PAS)
L = 100e-3 # Propagationssträcka (dvs avstånd mellan Plan 1 och 2)
lambda_noll = 633e-9 # Vakuumvåglängd för rött ljus från en HeNe-laser
n_medium = 1 # Brytningsindex för medium mellan Plan 1 och 2
k = 2*np.pi*n_medium/lambda_noll # K-vektorns längd
# Definera koordianter i plan 1
x = np.arange(-(N/2)*a, (N/2)*a, a) # Vektor med sampelpositioner i x-led
y = x # och y-led
X, Y = np.meshgrid(x, y) # Koordinatmatriser med x- och y-värdet i varje sampelposition
R = np.sqrt(X**2 + Y**2) # Avståndet till origo för varje sampelpunkt
# Definera lins och cirkulär aperatur
f_lins = 100e-3 # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins)) # Transmissionsfunktion för en lins (linsen är TOK)
D_aperture = 2e-3 # Diameter för apertur
T_aperture = R < (D_aperture/2) # Transmissionsfunktion för en cirkulär apertur ("pupill")
# Definera fält i plan 1
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'
omega1 = 1e-3 # 1/e2-radie (för intensiteten, dvs 1/e-radie för amplituden) för infallande Gaussiskt fält
E1_in_gauss = np.exp(-R**2/omega1**2) # Infallande fält: Gaussiskt med plana vågfronter och normalinfall (dvs konstant fas, här=0)
E1_in_konst = np.ones(X.shape) # Infallande fält: Konstant i hela plan 1. np.ones(X.shape) ger en matris fylld med ettor som har samma storlek som X
E1_gauss = E1_in_gauss*T_lins # Fältet i Plan 1 (precis efter linsen) för gaussisk stråle
E1_cirkular = E1_in_konst* T_lins* T_aperture # Fältet i Plan 1 (precis efter linsen) för konstant fält som passerat genom cirkulär apertur
E1 = E1_gauss # Välj fall!
I1 = np.abs(E1)**2 # Intensiteten är prop mot kvadraten på fältets amplitud (normalt struntar man i proportionalitetskonstanten)
I1_norm = I1/np.max(I1) # Då vi inte är intresserade av det absoluta värdet av intensiteten så kan det vara trevligt att noramlizera intensiten innan vi plottar
# Plotta intensitet i plan 1
x_mm = x*1e3
y_mm = y*1e3
plt.figure()
image = plt.imshow(I1_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 1. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 1
plt.figure()
image = plt.imshow(np.angle(E1), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Fas i plan 1. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Propagera till plan 2
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'
L = 1e-2
E2 = PAS(E1, L, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2) # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2 för L = 1 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Fas i plan 2 för L = 1 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
L = 5e-2
E2 = PAS(E1, L, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2) # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2 för L = 5 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Fas i plan 2 för L = 5 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
L = 8e-2
E2 = PAS(E1, L, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2) # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2 för L = 8 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Fas i plan 2 för L = 8 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
L = 10e-2
E2 = PAS(E1, L, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2) # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2 för L = 10 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Fas i plan 2 för L = 10 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# KOD
# Propagera till plan 2
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'
f_1 = 10e-2
f_2 = 1
def D_spot_gauss(f):
f_lins = f # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
E1 = E1_in_gauss*T_lins
E2 = PAS(E1, f, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2)
I2_norm_flat = I2_norm.flatten() #Vi flattenar matriserna så att vi enkelt kan hitte det R där I2_norm är närmast 1/e^2
R_flat = R.flatten()
index = np.argmin(np.abs(I2_norm_flat - 1/np.exp(2)))
omega = R_flat[index]
return 2*omega
C = D_spot_gauss(f_1) *2e-3 /(633e-9*f_1) # Vi löser ut C ur tumregelekvationen
print(C)
1.258644232225837
# Plotta intensitet i plan 2
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
C = D_spot_gauss(f_2) *2e-3 /(633e-9 *f_2)
print(C)
1.2835703001579646
f_lins = f_2 # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
E1 = E1_in_gauss*T_lins
E2 = PAS(E1, f_2, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2)
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
1.258644232225837 för f = 10 cm och 1.2835703001579646 för f = 1m så vi skulle säga att ett mer exakt värde skulle ligga vid 1.27
def D_spot_cirk(f):
f_lins = f # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
E1 = E1_in_konst*T_lins*T_aperture
E2 = PAS(E1, f, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2)
I2_norm_flat = I2_norm.flatten() #Vi flattenar matriserna så att vi enkelt kan hitte det R där I2_norm är närmast 1/e^2
R_flat = R.flatten()
index = np.argmin(np.abs(I2_norm_flat - 1/np.exp(2)))
omega = R_flat[index]
return 2*omega
C1 = D_spot_cirk(f_1) *2e-3 /(lambda_noll*f_1)
print(C1)
C2 = D_spot_cirk(f_2) *2e-3 /(lambda_noll*f_2)
print(C2)
1.6558560259922914 1.6272758416123274
1.6558560259922914 för f = 10cm och 1.6272758416123274 för f = 1m
f_lins = f_2 # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
E1 = E1_in_konst*T_lins*T_aperture
E2 = PAS(E1, f_2, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2)
plt.plot(Y[:, 0], I2_norm[:, I2_norm.shape[1]//2]) # Intensitetsprofil längs en linje
plt.xlabel("Y-position (m)")
plt.ylabel("Intensitet")
plt.title("Intensitet i fjärrfältet för f = 1m")
plt.grid(True)
plt.show()
f_lins = f_1 # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
E1 = E1_in_konst*T_lins*T_aperture
E2 = PAS(E1, f_1, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2)
plt.xlim(-0.0002, 0.0002)
plt.plot(Y[:, 0], I2_norm[:, I2_norm.shape[1]//2]) # Intensitetsprofil längs en linje
plt.xlabel("Y-position (m)")
plt.ylabel("Intensitet")
plt.title("Intensitet i fjärrfältet för f = 10 cm")
plt.grid(True)
plt.show()
graferna visar att intensitetsfördelningen är mycket brusigare för f = 1m och vi tror att detta kan vara ett tecken. Detta skulle även vara logiskt för alla fält som "åker" långa sträckor i ett homogent fält kan expandera vilket betyder att de kan expandera tillräckligt för att punkter ska hamna utanför det numeriska fönstret, alltså N b x N b.
T_lins2 = np.exp(-1j*k*R**2/(2*f_2))
E1_hermitgauss = E1_in_gauss*X
E2 = PAS(E1_hermitgauss, 10e-2, N, a, lambda_noll, n_medium) # Propagera med vår PAS funktion
I2 = np.abs(E2)**2 # Intesitet i plan 2
I2_norm = I2/np.max(I2)
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Fas i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
#Vi har kollat med många fokallängder mellan 0 och 1 meter och formen bibehålls.
# Ladda in DOE såhär!
DOE = np.flip(io.loadmat('T_DOE_gen2.mat'))["T_DOE_gen2"]
print(DOE)
[[0 0 0 ... 0 0 0] [0 0 0 ... 0 0 0] [0 0 0 ... 0 0 0] ... [0 0 0 ... 0 0 0] [0 0 0 ... 0 0 0] [0 0 0 ... 0 0 0]]
# Skapa det elektriska fältet som träffar ögat och propagera med pas
T_lins3 = np.exp(-1j*k*R**2/(2*0.02))
E1_in = E1_in_konst* T_lins3*DOE
E2 = PAS(E1_in, 0.02, N, a, 633e-9, 1)
# Fältet som har propagerats genom DOE och 'de vises lins' är smidigt att
# plotta i logaritmisk skala för att bilden ska vara tydlig.
# Alltså plotta fältet i plan 2 så här!
I2 = np.abs(E2)**2
I2_norm = np.log(I2/np.max(I2)) # Log av den normaliserade intensiteten i plan 2
x_mm = x*1e3
y_mm = y*1e3
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
<matplotlib.colorbar.Colorbar at 0x159c892e0d0>
T_lins_vise = np.exp(-1j*k*R**2/(2*0.15)) # Vi testade en del styrkor mellan 0 och 0.5, 0.15 var bäst!
E1_in = E1_in_konst* T_lins3*DOE*T_lins_vise
E2 = PAS(E1_in, 0.02, N, a, 633e-9, 1)
I2 = np.abs(E2)**2
I2_norm = np.log(I2/np.max(I2)) # Log av den normaliserade intensiteten i plan 2
x_mm = x*1e3
y_mm = y*1e3
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
<matplotlib.colorbar.Colorbar at 0x159c8294750>
Make Newton great again!
0.15 var den bästa styrkan och budskapet var <3 Fresnel